Wastewater Surviellence for SARS-CoV-2 in Athens, GA
Lipp Laboratory, Department of Environmental Health Science, University of Georgia
2020-10-30
Load all required libraries.
library(tidyverse)## Warning: package 'tidyverse' was built under R version 3.6.3## -- Attaching packages ---------------------------------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.0
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## Warning: package 'ggplot2' was built under R version 3.6.3## Warning: package 'tibble' was built under R version 3.6.3## Warning: package 'readr' was built under R version 3.6.3## Warning: package 'dplyr' was built under R version 3.6.3## Warning: package 'forcats' was built under R version 3.6.3## -- Conflicts ------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(plotly)## Warning: package 'plotly' was built under R version 3.6.3##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(broom)## Warning: package 'broom' was built under R version 3.6.3Read in raw data from RDS.
raw_data <- readRDS("./n1_n2_cleaned_cases.rds")Make a few small modifications to names and data for visualizations.
final_data <- raw_data %>% mutate(log_copy_per_L = log10(mean_copy_num_L)) %>%
rename(Facility = wrf) %>%
mutate(Facility = recode(Facility,
"NO" = "WRF A",
"MI" = "WRF B",
"CC" = "WRF C"))Seperate the data by gene target to ease layering in the final plot
#make three data layers
only_positives <<- subset(final_data, (!is.na(final_data$Facility)))
only_n1 <- subset(only_positives, target == "N1")
only_n2 <- subset(only_positives, target == "N2")
only_background <<-final_data %>%
select(c(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke, cases_per_100000_clarke)) %>%
group_by(date) %>% summarise_if(is.numeric, mean)
#specify fun colors
background_color <- "#7570B3"
seven_day_ave_color <- "#E6AB02"
marker_colors <- c("N1" = '#1B9E77',"N2" ='#D95F02')#remove facilty C for now
only_n1 <- only_n1[!(only_n1$Facility == "WRF C"),]
only_n2 <- only_n2[!(only_n2$Facility == "WRF C"),]Build the main plot
#first layer is the background epidemic curve
p1 <- only_background %>%
plotly::plot_ly() %>%
plotly::add_trace(x = ~date, y = ~new_cases_clarke,
type = "bar",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Daily Cases: ', new_cases_clarke),
alpha = 0.5,
name = "Daily Reported Cases",
color = background_color,
colors = background_color,
showlegend = FALSE) %>%
layout(yaxis = list(title = "Clarke County Daily Cases", showline=TRUE)) %>%
layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
#renders the main plot layer two as seven day moving average
p1 <- p1 %>% plotly::add_trace(x = ~date, y = ~X7_day_ave_clarke,
type = "scatter",
mode = "lines",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Seven-Day Moving Average: ', X7_day_ave_clarke),
name = "Seven Day Moving Average Athens",
line = list(color = seven_day_ave_color),
showlegend = FALSE)
#renders the main plot layer three as positive target hits
p2 <- plotly::plot_ly() %>%
plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Facility: ', Facility,
'</br> Target: ', target,
'</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
data = only_n1,
symbol = ~Facility,
marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Facility: ', Facility,
'</br> Target: ', target,
'</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
data = only_n2,
symbol = ~Facility,
marker = list(color = '#D95F02', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
layout(yaxis = list(title = "SARS CoV-2 Copies/L",
showline = TRUE,
type = "log",
dtick = 1,
automargin = TRUE)) %>%
layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
#adds the limit of detection dashed line
p2 <- p2 %>% plotly::add_segments(x = as.Date("2020-03-14"),
xend = ~max(date + 10),
y = 3571.429, yend = 3571.429,
opacity = 0.35,
line = list(color = "black", dash = "dash")) %>%
layout(annotations = list(x = as.Date("2020-03-28"), y = 3.8, xref = "x", yref = "y",
text = "Limit of Detection", showarrow = FALSE))
p1## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.## Warning: Ignoring 1 observations p2## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.p_combined <-
plotly::subplot(p2,p1, # plots to combine, top to bottom
nrows = 2,
heights = c(.6,.4), # relative heights of the two plots
shareX = TRUE, # plots will share an X axis
titleY = TRUE
) %>%
# create a vertical "spike line" to compare data across 2 plots
plotly::layout(
xaxis = list(
spikethickness = 1,
spikedash = "dot",
spikecolor = "black",
spikemode = "across+marker",
spikesnap = "cursor"
),
yaxis = list(spikethickness = 0)
)## Warning: Ignoring 1 observationsp_combinedSave the plot to pull into the index
save(p_combined, file = "./plotly_fig.rda")Save an htmlwidget for website embedding
htmlwidgets::saveWidget(p_combined, "plotly_fig.html")Build loess smoothing figures figures
#create smoothing data frames
#n1
smooth_n1 <- only_n1 %>% select(-c(Facility)) %>%
group_by(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke, cases_per_100000_clarke) %>%
summarize(sum_copy_num_L = sum(mean_total_copies)) %>%
ungroup() %>%
mutate(log_sum_copies_L = log10(sum_copy_num_L)) %>%
mutate(target = "N1")## `summarise()` regrouping output by 'date', 'cases_cum_clarke', 'new_cases_clarke', 'X7_day_ave_clarke' (override with `.groups` argument)#n2
smooth_n2 <- only_n2 %>% select(-c(Facility)) %>%
group_by(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke, cases_per_100000_clarke) %>%
summarize(sum_copy_num_L = sum(mean_total_copies)) %>%
ungroup() %>%
mutate(log_sum_copies_L = log10(sum_copy_num_L)) %>%
mutate(target = "N2")## `summarise()` regrouping output by 'date', 'cases_cum_clarke', 'new_cases_clarke', 'X7_day_ave_clarke' (override with `.groups` argument)#add trendlines
#extract data from geom_smooth
#n1 extract
# *********************************span 0.6***********************************
#*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_n1 <- ggplot(smooth_n1, aes(x = date, y = log_sum_copies_L)) +
stat_smooth(aes(outfit=fit_n1<<-..y..), method = "loess", color = '#1B9E77',
span = 0.6, n = 134)## Warning: Ignoring unknown aesthetics: outfit#n2 extract
extract_n2 <- ggplot(smooth_n2, aes(x = date, y = log_sum_copies_L)) +
stat_smooth(aes(outfit=fit_n2<<-..y..), method = "loess", color = '#1B9E77',
span = 0.6, n = 134)## Warning: Ignoring unknown aesthetics: outfit#look at the fits to align dates and total observations
#n1
extract_n1## `geom_smooth()` using formula 'y ~ x'
fit_n1## [1] 11.05203 11.21442 11.37213 11.52460 11.67132 11.81174 11.94531 12.07150
## [9] 12.19043 12.30282 12.40909 12.50966 12.60495 12.69536 12.78133 12.86125
## [17] 12.93358 12.99874 13.05719 13.10936 13.15570 13.19665 13.23265 13.26413
## [25] 13.29156 13.31535 13.33597 13.35384 13.36942 13.37309 13.35692 13.32380
## [33] 13.27667 13.21844 13.15201 13.08031 13.00626 12.93277 12.86275 12.79912
## [41] 12.74480 12.70271 12.67575 12.64827 12.60699 12.55845 12.50915 12.46563
## [49] 12.43440 12.42198 12.42463 12.43384 12.44896 12.46937 12.49445 12.52357
## [57] 12.55609 12.59138 12.62883 12.66779 12.70764 12.74775 12.78749 12.82624
## [65] 12.87378 12.93840 13.01719 13.10726 13.20570 13.30961 13.41609 13.52224
## [73] 13.62515 13.72192 13.80966 13.88545 13.94641 13.98961 14.02604 14.06646
## [81] 14.10767 14.14651 14.17977 14.20428 14.21684 14.21496 14.19949 14.17225
## [89] 14.13510 14.08988 14.03842 13.98257 13.92417 13.86506 13.80708 13.75208
## [97] 13.70188 13.65835 13.62331 13.59148 13.55679 13.51983 13.48124 13.44163
## [105] 13.40161 13.36180 13.32282 13.28529 13.24983 13.21704 13.18756 13.16199
## [113] 13.14095 13.12157 13.10116 13.08056 13.06066 13.04232 13.02641 13.01380
## [121] 13.00389 12.99553 12.98873 12.98348 12.97980 12.97768 12.97713 12.97815
## [129] 12.98075 12.98493 12.99070 12.99806 13.00701 13.01757#n2
extract_n2## `geom_smooth()` using formula 'y ~ x'
fit_n2## [1] 10.84581 11.04207 11.23287 11.41779 11.59640 11.76825 11.93292 12.08998
## [9] 12.23949 12.38202 12.51792 12.64752 12.77117 12.88920 13.00195 13.10821
## [17] 13.20676 13.29796 13.38212 13.45960 13.53073 13.59585 13.65530 13.70942
## [25] 13.75855 13.80302 13.84318 13.87936 13.91190 13.93284 13.93559 13.92256
## [33] 13.89619 13.85888 13.81305 13.76114 13.70554 13.64869 13.59301 13.54090
## [41] 13.49480 13.45711 13.43027 13.39572 13.33938 13.27036 13.19778 13.13073
## [49] 13.07835 13.04972 13.03754 13.02833 13.02201 13.01847 13.01762 13.01937
## [57] 13.02362 13.03027 13.03924 13.05042 13.06371 13.07904 13.09629 13.11537
## [65] 13.14377 13.18731 13.24343 13.30957 13.38317 13.46167 13.54252 13.62314
## [73] 13.70097 13.77347 13.83806 13.89219 13.93329 13.95881 13.98006 14.00754
## [81] 14.03746 14.06604 14.08949 14.10403 14.10585 14.09270 14.06611 14.02801
## [89] 13.98032 13.92498 13.86392 13.79907 13.73236 13.66571 13.60106 13.54034
## [97] 13.48547 13.43839 13.40103 13.36759 13.33146 13.29341 13.25417 13.21449
## [105] 13.17512 13.13679 13.10026 13.06627 13.03556 13.00889 12.98699 12.97061
## [113] 12.96049 12.95317 12.94539 12.93820 12.93268 12.92990 12.93093 12.93684
## [121] 12.94690 12.95967 12.97516 12.99336 13.01428 13.03792 13.06428 13.09335
## [129] 13.12515 13.15966 13.19690 13.23686 13.27954 13.32494#assign fits to a vector
n1_trend <- fit_n1
n2_trend <- fit_n2
#extract y min and max for each
limits_n1 <- ggplot_build(extract_n1)$data## `geom_smooth()` using formula 'y ~ x'limits_n1 <- as.data.frame(limits_n1)
n1_ymin <- limits_n1$ymin
n1_ymax <- limits_n1$ymax
limits_n2 <- ggplot_build(extract_n2)$data## `geom_smooth()` using formula 'y ~ x'limits_n2 <- as.data.frame(limits_n2)
n2_ymin <- limits_n2$ymin
n2_ymax <- limits_n2$ymax
#reassign dataframes (just to be safe)
work_n1 <- smooth_n1
work_n2 <- smooth_n2
#fill in missing dates to smooth fits
work_n1 <- work_n1 %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_n1 <- work_n1$date
work_n2 <- work_n2 %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_n2 <- work_n2$date
#create a new smooth dataframe to layer
smooth_frame_n1 <- data.frame(date_vec_n1, n1_trend, n1_ymin, n1_ymax)
smooth_frame_n2 <- data.frame(date_vec_n2, n2_trend, n2_ymin, n2_ymax)#make plotlys
#plot smooth frames
p3 <- plotly::plot_ly() %>%
plotly::add_lines(x = ~date_vec_n1, y = ~n1_trend,
data = smooth_frame_n1,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_n1,
'</br> Median Log Copies: ', round(n1_trend, digits = 2),
'</br> Target: N1'),
line = list(color = '#1B9E77', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
plotly::add_lines(x = ~date_vec_n2, y = ~n2_trend,
data = smooth_frame_n2,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_n2,
'</br> Median Log Copies: ', round(n2_trend, digits = 2),
'</br> Target: N2'),
line = list(color = '#D95F02', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
plotly::add_ribbons(x ~date_vec_n1, ymin = ~n1_ymin, ymax = ~n1_ymax,
showlegend = FALSE,
opacity = 0.25,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_n1, #leaving in case we want to change
'</br> Max Log Copies: ', round(n1_ymax, digits = 2),
'</br> Min Log Copies: ', round(n1_ymin, digits = 2),
'</br> Target: N1'),
name = "",
line = list(color = '#1B9E77')) %>%
plotly::add_ribbons(x ~date_vec_n2, ymin = ~n2_ymin, ymax = ~n2_ymax,
showlegend = FALSE,
opacity = 0.25,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_n2, #leaving in case we want to change
'</br> Max Log Copies: ', round(n2_ymax, digits = 2),
'</br> Min Log Copies: ', round(n2_ymin, digits = 2),
'</br> Target: N2'),
name = "",
line = list(color = '#D95F02')) %>%
layout(yaxis = list(title = "Total Log SARS CoV-2 Copies",
showline = TRUE,
automargin = TRUE)) %>%
layout(xaxis = list(title = "Date")) %>%
plotly::add_segments(x = as.Date("2020-06-24"),
xend = as.Date("2020-06-24"),
y = ~min(n1_ymin), yend = ~max(n1_ymax),
opacity = 0.35,
name = "Bars Repoen",
hoverinfo = "text",
text = "</br> Bars Reopen",
"</br> 2020-06-24",
showlegend = FALSE,
line = list(color = "black", dash = "dash")) %>%
plotly::add_segments(x = as.Date("2020-07-09"),
xend = as.Date("2020-07-09"),
y = ~min(n1_ymin), yend = ~max(n1_ymax),
opacity = 0.35,
name = "Mask Mandate",
hoverinfo = "text",
text = "</br> Mask Mandate",
"</br> 2020-07-09",
showlegend = FALSE,
line = list(color = "black", dash = "dash")) %>%
plotly::add_segments(x = as.Date("2020-08-20"),
xend = as.Date("2020-08-20"),
y = ~min(n1_ymin), yend = ~max(n1_ymax),
opacity = 0.35,
name = "</br> Classes Begin",
"</br> 2020-08-20",
hoverinfo = "text",
text = "Classes Begin",
showlegend = FALSE,
line = list(color = "black", dash = "dash")) %>%
plotly::add_segments(x = as.Date("2020-10-03"),
xend = as.Date("2020-10-03"),
y = ~min(n1_ymin), yend = ~max(n1_ymax),
opacity = 0.35,
name = "</br> First Home Football Game",
"</br> 2020-10-03",
hoverinfo = "text",
text = "First Home Football Game",
showlegend = FALSE,
line = list(color = "black", dash = "dash")) %>%
plotly::add_markers(x = ~date, y = ~log_sum_copies_L,
data = smooth_n1,
hoverinfo = "text",
showlegend = FALSE,
text = ~paste('</br> Date: ', date,
'</br> Observed Log Copies: ', round(log_sum_copies_L, digits = 2)),
marker = list(color = '#1B9E77', size = 6, opacity = 0.65)) %>%
plotly::add_markers(x = ~date, y = ~log_sum_copies_L,
data = smooth_n2,
hoverinfo = "text",
showlegend = FALSE,
text = ~paste('</br> Date: ', date,
'</br> Observed Log Copies: ', round(log_sum_copies_L, digits = 2)),
marker = list(color = '#D95F02', size = 6, opacity = 0.65))
p3Create final trend plot by stacking with epidemic curve
smooth_extract <-
plotly::subplot(p3,p1, # plots to combine, top to bottom
nrows = 2,
heights = c(.6,.4), # relative heights of the two plots
shareX = TRUE, # plots will share an X axis
titleY = TRUE
) %>%
# create a vertical "spike line" to compare data across 2 plots
plotly::layout(
xaxis = list(
spikethickness = 1,
spikedash = "dot",
spikecolor = "black",
spikemode = "across+marker",
spikesnap = "cursor"
),
yaxis = list(spikethickness = 0)
)## Warning: Ignoring 1 observationssmooth_extractsave(smooth_extract, file = "./smooth_extract.rda")